/*
 *  lib-smv: Matrix - Vector Library targeted to Structural Mechanics Problems
 *  Copyright (C) 2006-2008  St.Brcic, Lj.Zugic
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
/** @file LinEqs.cpp
 *  @brief Implementation of class LinEqs - solution of linear algebraic equations
 *
 *  LINEAR EQUATIONS  
 *  A * x = b (and/or A * X = B)
 *
 *  Direct and iterative solutions
 *  Crout's LU decompostion
 *  Cholesky L LT decomposition
 *  Conjugate gradients (not yet)
 */

#include "LinEqs.h"
namespace smv_lib
{

//-------------------------------- constructors
//! empty (default) constructor
LinEqs::LinEqs(void)
{
  isMatrix_A  = false;
  isVector_b  = false;
  isMatrix_B  = false;
  isMore      = false;
  isNameGiven = false;
  isDirect    = false;
  isCrout     = false;
  isCholesky  = false;
  isSolved    = false;
}

//! constructor with a given matrix A and vector b - to solve equation A*x=b
LinEqs::LinEqs(const Matrix& AA, const Vector& bb, EQ_GROUP group)
    : A(AA), b(bb), EQ_Group(group)
{
  N = A.getRows();
  assert(N == A.getCols());
  isMatrix_A  = true;
  Name = A.getName();
  if (Name == "MATRIX")
    isNameGiven = false;
  else
    isNameGiven = true;
  A.Classify();
  if (EQ_Group == DIRECT) {
    switch (A.get_MatType()) {
        case DIAGONAL:
        case SYMMETRIC:
        case TRI_SYMM:
        EQ_Direct  = CHOLESKY;
        isCholesky = true;
        isCrout    = false;
        break;
        case NOT_CLASSIFIED:
        case GEN_TYPE:
        case TRIDIAG:
        EQ_Direct  = CROUT;
        isCrout    = true;
        isCholesky = false;
        break;
    }
    isDirect = true;
  }
  //--- right-hand side vector b
  assert (N == b.getLen());
  isVector_b = true;
  isMatrix_B = false;
  //--- solution vector x
  Vector* tmp = new Vector(N);
  assert(tmp != NULL);
  x = (*tmp);
  delete tmp;
  isSolved  = false;
  isMore    = false;
}

//! constructor with a given matrix A and a matrix B - solution of equations A*X=B
LinEqs::LinEqs(const Matrix& AA, const Matrix& BB, EQ_GROUP group)
    : A(AA), B(BB), EQ_Group(group)
{
  //--- matrix A (coefficient matrix)
  N = A.getRows();
  assert(N == A.getCols());
  isMatrix_A = true;
  Name = A.getName();
  if (Name == "MATRIX")
    isNameGiven = false;
  else
    isNameGiven = true;
  A.Classify();
  if (EQ_Group == DIRECT) {
    switch (A.get_MatType()) {
        case DIAGONAL:
        case SYMMETRIC:
        case TRI_SYMM:
        EQ_Direct  = CHOLESKY;
        isCholesky = true;
        isCrout    = false;
        break;
        case NOT_CLASSIFIED:
        case GEN_TYPE:
        case TRIDIAG:
        EQ_Direct  = CROUT;
        isCrout    = true;
        isCholesky = false;
        break;
    }
    isDirect = true;
  }     //--- matrix B (right-hand side)
  assert (N == B.getRows());
  M = B.getCols();
  isMatrix_B = true;
  isVector_b = false;
  //--- matrix X (solution of equations A * X = B)
  Matrix* tmp = new Matrix(N,M);
  assert(tmp != NULL);
  X = (*tmp);
  delete tmp;
  isSolved  = false;
  isMore    = true;
}

//! destructor
LinEqs::~LinEqs(void)
{
  if (isDirect == true) {
    if (isCrout == true)
      delete m_Crout;
    else if (isCholesky == true)
      delete m_Chol;
  }
}

//! set the coefficient matrix A (in the case of empty constructor)
void LinEqs::setMatrix_A(const Matrix& AA)
{
  //--- matrix A
  A = AA;
  N = A.getRows();
  assert(N == A.getCols());
  isMatrix_A  = true;
  Name = A.getName();
  if (Name == "MATRIX")
    isNameGiven = false;
  else
    isNameGiven = true;
  A.Classify();
  isVector_b    = false;
  isMatrix_B    = false;
  isSolved      = false;
  isMore        = false;
}

//! set the right-hand side vector b (Matrix A must be given first)
void LinEqs::setVector_b(const Vector& bb)
{
  if (isMatrix_A == false) {
    cerr << "\nGive the Matrix A first!";
    return;
  }
  b = bb;
  assert(N == b.getLen());
  isVector_b = true;
  isMatrix_B = false;
  //--- solution vector x
  Vector* tmp = new Vector(N);
  assert(tmp != NULL);
  x = (*tmp);
  delete tmp;
  isSolved      = false;
  isMore        = false;
}

//! set the right-hand side matrix B (Matrix A must be given first)
void LinEqs::setMatrix_B(const Matrix& BB)
{
  if (isMatrix_A == false) {
    cerr << "\nGive the Matrix A first!";
    return;
  }
  //--- matrix B (right-hand side)
  B = BB;
  assert (N == B.getRows());
  M = B.getCols();
  isMatrix_B = true;
  isVector_b = false;
  //--- matrix X (solution of equations A * X = B)
  Matrix* tmp = new Matrix(N,M);
  assert(tmp != NULL);
  X = (*tmp);
  delete tmp;
  isSolved      = false;
  isMore        = true;
}

//! set the name of the problem
void LinEqs::setName(const string name)
{
  Name = name;
  isNameGiven = true;
}

//! set the method of solution: default is direct method
void LinEqs::setMethod(EQ_GROUP group)
{
  EQ_Group = group;
  if (EQ_Group == DIRECT)
    isDirect = true;
  if (EQ_Group == ITERATIVE)
    isDirect = false;
}

//! set direct method of solution: default is Crout's LU
void LinEqs::setDirect(EQ_DIRECT direct)
{
  EQ_Direct = direct;
  if (EQ_Direct == CROUT) {
    isCrout    = true;
    isCholesky = false;
  }
  else if (EQ_Direct == CHOLESKY) {
    isCrout    = false;
    isCholesky = true;
  }
}

//================================= public working functions
//
//! solve the problem - find the solution of linear equations
void LinEqs::solveProblem()
{
  //--- coefficient matrix and the right-hand side must be given
  if (isMatrix_A != true) {
    cerr << "Coefficient matrix A is NOT given!\n";
    return;
  }
  if (isVector_b == true || isMatrix_B == true) {
    switch (EQ_Group) {
        case DIRECT:
           direct();
           break;
        case ITERATIVE:
           iterative();
           break;
    }
    isSolved = true;
    return;
  }
  else {
    cerr << "Right-hand side is NOT given!\n";
    return;
  }
}

//! get solution of linear equations A * x = b
Vector LinEqs::getSolution()
{
  assert(isSolved == true);
  return x;
}

//! get one solutions of the system of linear equations A * X = B
Vector LinEqs::getSolutionNo(int i)
{
  assert(isSolved == true);
  assert(i >=0 && i<M);
  return (X.getColData(i));
}

//! get all solutions of the system of linear equations A * X = B
Matrix LinEqs::getAllSolutions()
{
  assert(isSolved == true);
  return X;
}

//! write results to a file (with predefined name "LinEqs.txt")
void LinEqs::Write()
{
  string FileName;
  if (isNameGiven == false)
    FileName = "LinEqs.txt";
  else
    FileName = Name + ".txt";
  this->WriteTo(FileName.c_str());
}

//! write results to a given file
void LinEqs::WriteTo(const char* file)
{
  //--- open the file
  out.open(file,ios::out);
  assert(out != NULL);
  out << "===================== SOLUTION OF LINEAR EQUATIONS ========================" << endl << endl;
  out << "(1) INPUT MATRIX DATA: " << endl;
  out << "    - Name is:  " << Name << endl;
  out << "    - Order N = " << N << endl << endl;
  out << A << endl;
  out << giveMethod() << endl;
  if (isDirect == true) {
    out << giveDirectMethod() << endl;
  } else {
    out << giveIterative() << endl;
  }
  if (isSolved == true && isMore == false) {
    out << "(2) SOLUTION OF LINEAR EQUATIONS: " << endl << endl;
    out << "    Right-hand side = " << endl;
    out << b << endl;
    out << "    The solution is: " << endl;
    out << x << endl;
  }
  if (isSolved == true && isMore == true) {
    out << "(2) SOLUTION OF THE SYSTEM OF LINEAR EQUATIONS: " << endl << endl;
    for (int i=0; i<M; i++) {
      out << "   System number " << (i+1) << endl;
      out << "   Right-hand side = " << endl;
      out << B.getColData(i) << endl;
      out << "   The solution is: " << endl;
      out << X.getColData(i) << endl << endl;
    }
  }
  out.close();
}

//============================================ private functions
//! direct solution of linear equations
void LinEqs::direct()
{
  switch (EQ_Direct) {
      case CROUT:
         crout();
         break;
      case CHOLESKY:
         cholesky();
         break;
  }
}

//! Crout's LU decomp[osition (general matrix)
void LinEqs::crout()
{
  //--- call Crout (for general matrix)
  if (isVector_b == true)  {
    m_Crout = new Crout(A,b);
    assert (m_Crout != NULL);
  }
  else if (isMatrix_B == true) {
    m_Crout = new Crout(A,B);
    assert (m_Crout != NULL);
  }
  //--- decompose and obtain solution
  m_Crout->Decomposition();
  m_Crout->Solution();
  //--- transfer the solution to LinEqs class
  if (isVector_b == true)
    x = m_Crout->getSolution();
  else if (isMatrix_B == true)
    X = m_Crout->getAllSolutions();
}

//! Cholesky decomposition (symmetric matrix)
void LinEqs::cholesky()
{
  //--- matrix A MUST be symmetric for Cholesky
  //--- call Cholesky
  if (isVector_b == true)  {
    m_Chol = new Cholesky(A,b);           // standard decomposition A =L*LT
    assert (m_Chol != NULL);
  }
  else if (isMatrix_B == true) {
    m_Chol = new Cholesky(A,B);           // standard decomposition A =L*LT
    assert (m_Chol != NULL);
  }
  //--- decompose and obtain solution
  m_Chol->Decomposition();
  m_Chol->Solution();
  //--- transfer the solution to LinEqs class
  if (isVector_b == true)
    x = m_Chol->getSolution();
  else if (isMatrix_B == true)
    X = m_Chol->getAllSolutions();
}

//! give the method of solution as a string
const string LinEqs::giveMethod() const
{
  string info = "Method of solutiontion is ";
  switch (isDirect) {
      case true:
         info = info + "DIRECT\n";
         break;
      case false:
         info = info + "ITERATIVE\n";
         break;
  }
  return (info);
}

//! give the method of direct solution as a string
string LinEqs::giveDirectMethod()
{
  string info = "DIRECT SOLUTION is ";
  switch (isCrout) {
      case true:
         info = info + "CROUT'S LU DECOMPOSITION\n";
         break;
      case false:
         info = info + "CHOLESKY DECOMPOSITION\n";
         break;
  }
  return (info);
}

//! iterative solutions of linear equations
void LinEqs::iterative()
{
//  waits for better times!
}

//! give the method of iterative solution as a string
string LinEqs::giveIterative()
{
  string info = " Method of ITERATIVE SOLUTION is ";
  info = info + "CONJUGATE GRADIENTS\n";

  return (info);
}

}  // end of namespace smv_lib
//========================================= end of file
